#include "algorithm/analytical_model/algorithm/dnn/params/UNIZAR2.h"
#include <bit>
#define SOLAR_USE_CUDA_INTER_OP_DATA

#include <algorithm/analytical_model/CommonPara.h>
#include <algorithm/analytical_model/Constants.h>
#include <algorithm/analytical_model/HeliostatPara.h>
#include <algorithm/analytical_model/ScenePara.h>
#include <algorithm/analytical_model/SimulationPara.h>
#include <algorithm/analytical_model/algorithm/AllModels.h>
#include <algorithm/analytical_model/algorithm/ModelPara.h>
#include <algorithm/analytical_model/algorithm/cpu/VectorConvert.h>
#include <algorithm/analytical_model/algorithm/cuda/Utils.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/ConvolutionMethod.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/Data.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/ShadowBlockCalculator.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/SimplifiedMethod.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/grid/RectangleGrid.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/heliostat/RectangleHeliostat.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/receiver/CylinderReceiver.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/receiver/RectangleReceiver.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/utils/convert/AllConverts.h>
#include <algorithm/analytical_model/algorithm/dnn/params/AllParams.h>
#include <algorithm/analytical_model/algorithm/geometry/Utils.h>
#include <algorithm/analytical_model/algorithm/pipeline/AnalyticalPipeline.h>
#include <algorithm/analytical_model/algorithm/pipeline/AnalyticalPipelineParam.h>
#include <algorithm/analytical_model/algorithm/pipeline/Utils.h>
#include <algorithm/analytical_model/convolution/AllParas.h>
#include <algorithm/analytical_model/io/csv/WriteFluxMap.h>
#include <algorithm/analytical_model/io/json/converter/AllConverters.h>
#include <algorithm/analytical_model/io/json/transform/TransformJSON.h>
#include <algorithm/analytical_model/simplified/AllParas.h>
#include <algorithm/shared/data/scene/Heliostat.h>
#include <algorithm/shared/data/scene/Receiver.h>
#include <algorithm/shared/data/scene/Sun.h>
#include <algorithm/shared/data/topology/geometry/Cylinder3D.h>
#include <algorithm/shared/data/topology/geometry/Rectangle3D.h>
#include <algorithm/shared/tool/stringHash/StringHash.hpp>
#include <algorithm/shared/utils/Vec3.h>
#include <utils/Format.h>
#include <utils/ParaUtils.h>
#include <utils/Print.h>
#include <utils/Utils.h>
#include <utils/config.h>

#include <boost/asio.hpp>
#include <boost/iostreams/device/mapped_file.hpp>
#include <boost/iostreams/stream.hpp>
#include <boost/lockfree/queue.hpp>
#include <boost/mp11.hpp>
#include <boost/timer/timer.hpp>

#include <oneapi/tbb.h>

#include <cuda_runtime_api.h>

#include <array>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <memory>
#include <sstream>
#include <variant>
#include <vector>

// xxx_row_col: x（横坐标像素个数） y（纵坐标像素个数）
// 注意：线程池不能为static变量（有奇怪的bug，会跳过任务执行）

using namespace solar;

namespace solar
{
    using SimplifiedModelPara = std::variant<HFLCALPara, IHFLCALPara, BHFLCALPara, NEGPara>;
    using ConvolutionModelPara = std::variant<UNIZARPara, HuangPara, ICauchyPara, UNIZARPara2>;
    using DnnBasedModelPara = std::variant<NEGPara, UNIZARPara2>;
    using SimplifiedModelParaTypeList =
        boost::mp11::mp_list<HFLCALPara, IHFLCALPara, BHFLCALPara, NEGPara>;
    using ConvolutionModelParaTypeList =
        boost::mp11::mp_list<UNIZARPara, HuangPara, ICauchyPara, UNIZARPara2>;
    using DnnBasedModelParaTypeList = boost::mp11::mp_list<NEGPara, UNIZARPara2>;
    using DnnBasedModelInputTypeList = boost::mp11::mp_list<dnn::NEGInput, dnn::UNIZAR2Input>;
    using DnnBasedModelOutputTypeList = boost::mp11::mp_list<dnn::NEGOutput, dnn::UNIZAR2Output>;
    // namespace internal
    // {
    //     static boost::asio::thread_pool preprocessing_thread_pool(3);
    // }
} // namespace solar

namespace solar::internal
{
    static auto getHeliostatRowCol(Geometry* geometry, const Vec2& heliostat_pixel_length) -> Int2
    {
        // std::cout << solar::format("heliostat_pixel_length: {}, {}\n",
        // heliostat_pixel_length.x(), heliostat_pixel_length.y());
        auto geometry_heliostat_type = geometry->getType();
        switch (geometry_heliostat_type)
        {
        case GeometryType::kGeometryRectangle:
        {
            auto* rectangle_geometry = dynamic_cast<Rectangle*>(geometry);
            auto size = SOLAR_GET_VALUE(Rectangle, size_, *rectangle_geometry);
            auto col_count = size.x() / heliostat_pixel_length.x();
            auto rol_count = col_count;
            if (solar::hasValue(heliostat_pixel_length.y()))
            {
                rol_count = size.y() / heliostat_pixel_length.y();
            }
            return Int2(col_count, rol_count);
        }
        case GeometryType::kGeometryRectangle3D:
        {
            auto* rectangle3d_geometry = dynamic_cast<Rectangle3D*>(geometry);
            auto size = (SOLAR_GET_VALUE(Rectangle3D, size_, *rectangle3d_geometry)).data();
            auto col_count = size[0] / heliostat_pixel_length.x();
            auto rol_count = col_count;
            if (solar::hasValue(heliostat_pixel_length.y()))
            {
                rol_count = size[2] / heliostat_pixel_length.y();
            }
            return Int2(col_count, rol_count);
        }
        default:
        {
            const auto error_string =
                ::solar::format("ERROR: In function {}:\n {} is not implemented",
                                PRETTY_FUNCTION_NAME, solar::getEnumName(geometry_heliostat_type));
            solar::println(error_string);
            throw error_string;
        }
        }
    }
} // namespace solar::internal

auto solar::runAnalyticalModelPipeline(AnalyticalPipelineParam* params, cudaStream_t stream)
    -> std::tuple<float*, cuda::inter_op::Int2>
{
    boost::asio::thread_pool preprocessing_thread_pool(3);

    auto receiver_pixel_size =
        SOLAR_GET_VALUE(SimulationPara, receiver_pixel_size_, params->simulation_para_);
    auto micro_heliostat_pixel_size =
        SOLAR_GET_VALUE(SimulationPara, micro_heliostat_pixel_size_, params->simulation_para_);

    // Assume there's only one receiver

    auto receiver_geometry = SOLAR_GET_VALUE(Receiver, geometry_, params->receivers_[0]);
    // auto receiver_center = solar::getCenter(receiver_geometry.get());
    // JSON中缺少相关信息（经纬度、时间等），暂时这么写
    // 计算代码见 SceneRendererAdapter::increaseTime()
    auto sun_position = SOLAR_GET_VALUE(Sun, sun_position_, params->sun_);
    const auto azimuth_degree = SOLAR_GET_VALUE(SunPosition, azimuth_, *sun_position);
    const auto elevation_degree = SOLAR_GET_VALUE(SunPosition, elevation_, *sun_position);
    Float3 sunray_dir = cpu::toFloat3(-angleToDirection(azimuth_degree, elevation_degree));

    // Configure grid

    cuda::inter_op::RectangleGrid grid;
    solar::cuda::inter_op::convertGrid(std::addressof(params->rectangle_grid_),
                                       std::addressof(grid));
    grid.setNumHeliostats(params->heliostats_.size());
    cuda::inter_op::ShadowBlockCalculator shadow_block_calculator;
    auto heliostat_tmp_geometry = SOLAR_GET_VALUE(Heliostat, geometry_, params->heliostats_[0]);
    const auto helio_row_col =
        internal::getHeliostatRowCol(heliostat_tmp_geometry.get(), micro_heliostat_pixel_size);
    // std::cout << solar::format("heliostat row_col: {}, {}\n", helio_row_col.x, helio_row_col.y);
    shadow_block_calculator.setHelioRowCol(helio_row_col);
    shadow_block_calculator.setSunRayDirection(sunray_dir);

    // Rectangle Heliostat only for now.
    // Although RayTracing Pipeline supports other types of heliostats, the only supported heliostat
    // type in analytical method is rectangle heliostat for now.
    // Maybe you can use multiple rectangle
    // heliostats to simulate a focused heliostat in analytical method [Not tested]

    std::vector<cuda::inter_op::RectangleHeliostat> rectangle_heliostats;
    std::vector<Float3> all_heliostat_vertices;
    std::vector<Float3> heliostat_positions(params->heliostats_.size());
    std::vector<Float3> heliostat_image_plane_axis;
    if (params->simulation_option_ == SimulationOption::Whole)
    {
        if (params->model_type_ == ModelType::kSimplified)
        {
            // 定日镜成像平面 x/z 轴
            heliostat_image_plane_axis.resize(params->heliostats_.size() * 2);
        }
        if (params->model_type_ == ModelType::kConvolution)
        {
            // 定日镜成像平面 x/y/z 轴
            heliostat_image_plane_axis.resize(params->heliostats_.size() * 3);
        }
    }
    rectangle_heliostats.resize(params->heliostats_.size());
    all_heliostat_vertices.resize(params->heliostats_.size() * 4);
    SOLAR_OMP_PARALLEL_FOR_SIMD
    for (int i = 0; i < params->heliostats_.size(); i++)
    {
        solar::cuda::inter_op::convertHeliostat(std::addressof(params->heliostats_[i]),
                                                std::addressof(rectangle_heliostats[i]));
        rectangle_heliostats[i].setRowAndColumn({1, 1});
        rectangle_heliostats[i].setBelongingGridId(0);
        auto heliostat_geometry = SOLAR_GET_VALUE(Heliostat, geometry_, params->heliostats_[i]);
        auto focus_center = getFocusCenter(heliostat_geometry.get(), receiver_geometry.get());
        rectangle_heliostats[i].setNormalAndRotate(cpu::toFloat3(focus_center), sunray_dir);
        auto j = i * 4;
        const auto current_vertices = rectangle_heliostats[i].getVertices();
        std::copy(std::addressof(current_vertices[0]), std::addressof(current_vertices[4]),
                  std::addressof(all_heliostat_vertices[j]));
    }
    Float3* d_all_heliostat_vertices = nullptr;
    if (params->simulation_option_ == SimulationOption::Whole)
    {
        if (params->model_type_ == ModelType::kConvolution)
        {
            solar::cuda::cpu2gpuAsync(d_all_heliostat_vertices, all_heliostat_vertices.data(),
                                      all_heliostat_vertices.size(),
                                      shadow_block_calculator.getCudaStream());
        }
    }

    // 高斯模型：x/z轴 卷积模型：xyz轴
    Float3* d_image_plane_axes = nullptr;
    if (params->simulation_option_ == SimulationOption::Whole)
    {
        CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(d_image_plane_axes),
                                          sizeof(Float3) * heliostat_image_plane_axis.size(),
                                          shadow_block_calculator.getCudaStream()));
    }

    boost::asio::post(preprocessing_thread_pool,
                      [&]()
                      {
                          SOLAR_OMP_PARALLEL_FOR_SIMD
                          for (int i = 0; i < params->heliostats_.size(); i++)
                          {
                              heliostat_positions[i] = rectangle_heliostats[i].getPosition();
                              if (params->simulation_option_ == SimulationOption::Whole)
                              {
                                  if (params->model_type_ == ModelType::kSimplified)
                                  {
                                      // 定日镜成像平面的x/z轴
                                      std::tie(heliostat_image_plane_axis[i * 2], std::ignore,
                                               heliostat_image_plane_axis[i * 2 + 1]) =
                                          rectangle_heliostats[i].getImagePlaneAxis();
                                  }
                                  if (params->model_type_ == ModelType::kConvolution)
                                  {
                                      std::tie(heliostat_image_plane_axis[i * 3],
                                               heliostat_image_plane_axis[i * 3 + 1],
                                               heliostat_image_plane_axis[i * 3 + 2]) =
                                          rectangle_heliostats[i].getImagePlaneAxis();
                                  }
                              }
                          }
                          if (params->simulation_option_ == SimulationOption::Whole)
                          {
                              solar::cuda::cpu2gpuAsync(d_image_plane_axes,
                                                        heliostat_image_plane_axis.data(),
                                                        heliostat_image_plane_axis.size(),
                                                        shadow_block_calculator.getCudaStream());
                          }
                      });

    std::vector<SimplifiedModelPara> simplified_model_paras;
    std::vector<ConvolutionModelPara> convolution_model_paras;
    if (params->model_type_ == ModelType::kSimplified)
    {
        simplified_model_paras.resize(params->heliostats_.size());
    }
    if (params->model_type_ == ModelType::kConvolution)
    {
        convolution_model_paras.resize(params->heliostats_.size());
    }

    cudaStream_t simulation_stream = stream;
    if (simulation_stream == nullptr)
    {
        CHECK_CUDA_ERRORS(cudaStreamCreate(std::addressof(simulation_stream)));
    }
    float* d_dnn_outputs = nullptr;
    bool use_dnn_based_model = false;

    static const auto compute_sigma = [&]<typename T>(const T*)
    {
        static_assert(!boost::mp11::mp_contains<DnnBasedModelPara, T>::value);
        SOLAR_OMP_PARALLEL_FOR_SIMD
        for (int i = 0; i < rectangle_heliostats.size(); i++)
        {
            convertHeliostat(std::addressof(rectangle_heliostats[i]),
                             std::addressof(params->heliostats_[i]));
            T current_model_para;
            solar::computeSigma(
                std::addressof(params->simulation_para_), std::addressof(params->scene_para_),
                std::addressof(params->heliostat_paras_[i]), std::addressof(params->common_para_),
                std::addressof(params->sun_), std::addressof(params->heliostats_[i]),
                std::addressof(params->receivers_[0]), std::addressof(current_model_para));
            if constexpr (boost::mp11::mp_contains<SimplifiedModelPara, T>::value)
            {
                simplified_model_paras[i] = current_model_para;
            }
            if constexpr (boost::mp11::mp_contains<ConvolutionModelPara, T>::value)
            {
                convolution_model_paras[i] = current_model_para;
            }
        }
    };

    static const auto compute_sigma_dnn = [&]<typename T>(const T*)
    {
        static_assert(boost::mp11::mp_contains<DnnBasedModelPara, T>::value);
        use_dnn_based_model = true;
        constexpr auto type_index = boost::mp11::mp_find<DnnBasedModelParaTypeList, T>::value;
        using model_output_type = boost::mp11::mp_at_c<DnnBasedModelOutputTypeList, type_index>;
        const auto model_output_size = sizeof(model_output_type) * params->heliostats_.size();
        CHECK_CUDA_ERRORS(
            cudaMallocAsync(std::addressof(d_dnn_outputs), model_output_size, simulation_stream));
        const auto compute_sigma_callback = [&](void* ptr)
        {
            CHECK_CUDA_ERRORS(cudaMemcpyAsync(d_dnn_outputs, ptr, model_output_size,
                                              cudaMemcpyHostToDevice, simulation_stream));
        };
        FuncConverter<void(void*)>::func_ = compute_sigma_callback;
        ComputeSigmaCallback compute_sigma_callback_func = FuncConverter<void(void*)>::callback;
        SOLAR_OMP_PARALLEL_FOR_SIMD
        for (int i = 0; i < rectangle_heliostats.size(); i++)
        {
            convertHeliostat(std::addressof(rectangle_heliostats[i]),
                             std::addressof(params->heliostats_[i]));
        }
        std::vector<T> model_paras(params->heliostats_.size());
        solar::computeSigma(std::addressof(params->simulation_para_),
                            std::addressof(params->scene_para_), params->heliostat_paras_.data(),
                            std::addressof(params->common_para_), std::addressof(params->sun_),
                            params->heliostats_.data(), std::addressof(params->receivers_[0]),
                            params->heliostats_.size(), model_paras.data(),
                            compute_sigma_callback_func);
        SOLAR_OMP_PARALLEL_FOR_SIMD
        for (int i = 0; i < model_paras.size(); i++)
        {
            const auto current_model_para = model_paras[i];
            if constexpr (boost::mp11::mp_contains<SimplifiedModelPara, T>::value)
            {
                simplified_model_paras[i] = current_model_para;
            }
            if constexpr (boost::mp11::mp_contains<ConvolutionModelPara, T>::value)
            {
                convolution_model_paras[i] = current_model_para;
            }
        }
    };

    boost::asio::post(preprocessing_thread_pool,
                      [&]()
                      {
                          // demo only
                          // JSON中暂时没有这些数据
                          SOLAR_OMP_PARALLEL_FOR_SIMD
                          for (int i = 0; i < params->heliostat_paras_.size(); i++)
                          {
                              SOLAR_SET_VALUE(HeliostatPara, sigma_tracking_, Vec2(0.0, 0.0),
                                              params->heliostat_paras_[i])
                              //   SOLAR_SET_VALUE(HeliostatPara, sigma_slope_error_, Vec2(0.001,
                              //   0.001),
                              //                   params->heliostat_paras_[i])
                          }

                          if (params->model_type_ == ModelType::kSimplified)
                          {
                              switch (std::get<SimplifiedModelType>(params->detailed_model_type_))
                              {
                              case SimplifiedModelType::kHFLCAL:
                              {
                                  HFLCALPara* ptr = nullptr;
                                  compute_sigma(ptr);
                                  break;
                              }
                              case SimplifiedModelType::kiHFLCAL:
                              {
                                  IHFLCALPara* ptr = nullptr;
                                  compute_sigma(ptr);
                                  break;
                              }
                              case SimplifiedModelType::kbHFLCAL:
                              {
                                  BHFLCALPara* ptr = nullptr;
                                  compute_sigma(ptr);
                                  break;
                              }
                              case SimplifiedModelType::kNEG:
                              {
                                  NEGPara* ptr = nullptr;
                                  compute_sigma_dnn(ptr);
                                  break;
                              }
                              }
                          }
                          if (params->model_type_ == ModelType::kConvolution)
                          {
                              switch (std::get<ConvolutionModelType>(params->detailed_model_type_))
                              {
                              case ConvolutionModelType::kUNIZAR:
                              {
                                  UNIZARPara* ptr = nullptr;
                                  compute_sigma(ptr);
                                  break;
                              }
                              case ConvolutionModelType::kHuang:
                              {
                                  HuangPara* ptr = nullptr;
                                  compute_sigma(ptr);
                                  break;
                              }
                              case ConvolutionModelType::kiCauchy:
                              {
                                  ICauchyPara* ptr = nullptr;
                                  compute_sigma(ptr);
                                  break;
                              }
                              case ConvolutionModelType::kUNIZAR2:
                              {
                                  UNIZARPara2* ptr = nullptr;
                                  compute_sigma_dnn(ptr);
                                  break;
                              }
                              }
                          }
                      });

    grid.setBelongingReceiverIndex(0);
    grid.setStartHeliostatIndex(0);
    grid.init();
    grid.gridHelioMatchGPU(rectangle_heliostats);
    solar::println("grid.gridHelioMatchGPU finished");
    // boost::asio::post(preprocessing_thread_pool,
    //                   [&]() { grid.gridHelioMatchGPU(rectangle_heliostats); });
    shadow_block_calculator.setGrid(&grid);
    shadow_block_calculator.init(all_heliostat_vertices);
    // boost::asio::post(preprocessing_thread_pool,
    //                   [&]() { shadow_block_calculator.init(all_heliostat_vertices); });
    // shadow_block_calculator.init(all_heliostat_vertices);
    // preprocessing_thread_pool.join();
    // cudaStreamSynchronize(grid.getCudaStream());
    solar::println("before shadow_block_calculator.calculate");
    shadow_block_calculator.calculate();

    auto h_shadow_block_center_bias_factor =
        std::make_shared<Float4[]>(rectangle_heliostats.size());
    auto* h_shadow_block_center_bias_factor_ptr = h_shadow_block_center_bias_factor.get();

    Float3* d_heliostat_positions = nullptr;
    auto* h_heliostat_positions = heliostat_positions.data();
    solar::cuda::cpu2gpuAsync(d_heliostat_positions, h_heliostat_positions,
                              rectangle_heliostats.size(), shadow_block_calculator.getCudaStream());

    if (params->model_type_ == ModelType::kSimplified)
    {

        // 定日镜阴影遮挡结果
        // x,y,z 为重心偏移后的中心坐标
        // w 为阴影遮挡比例
        shadow_block_calculator.calculateCenterBias(d_heliostat_positions);
    }

    static constexpr char sigma_t[] = "sigma_";
    static constexpr char sigma_x_t[] = "sigma_x_";
    static constexpr char sigma_y_t[] = "sigma_y_";
    static constexpr char x3_t[] = "x3_";
    static constexpr char x0_t[] = "x0_";
    static constexpr char y1_t[] = "y1_";
    static constexpr char y0_t[] = "y0_";

    static constexpr char diffusion_coeff_t[] = "diffusion_coeff_";
    static constexpr char flux_param_t[] = "flux_param_";
    static const auto get_single_sigma_params =
        [&]<typename T>(T* para) -> std::tuple<double, double, double>
    {
        const auto sigma = para->*solar::getParaPtr<T, sigma_t>();
        const auto flux_param = para->*solar::getParaPtr<T, flux_param_t>();
        return {sigma, sigma, flux_param};
    };
    static const auto get_double_sigma_params =
        [&]<typename T>(T* para) -> std::tuple<double, double, double>
    {
        const auto sigma_x = para->*solar::getParaPtr<T, sigma_x_t>();
        const auto sigma_y = para->*solar::getParaPtr<T, sigma_y_t>();
        const auto flux_param = para->*solar::getParaPtr<T, flux_param_t>();
        return {sigma_x, sigma_y, flux_param};
    };
    static const auto get_convolution_area_rectangle =
        [&]<typename T>(T* para) -> std::tuple<double, double, double, double>
    {
        const auto x3 = para->*solar::getParaPtr<T, x3_t>();
        const auto x0 = para->*solar::getParaPtr<T, x0_t>();
        const auto y1 = para->*solar::getParaPtr<T, y1_t>();
        const auto y0 = para->*solar::getParaPtr<T, y0_t>();
        // 右上 左上 左下 左上
        return {x3, x0, y1, y0};
    };

    // returns (sigma_x, sigma_y, flux_param)
    static const auto get_sigmas = [&](int index) -> std::tuple<double, double, double>
    {
        if (params->model_type_ == ModelType::kSimplified)
        {
            if (auto* hflcal_para =
                    std::get_if<HFLCALPara>(std::addressof(simplified_model_paras[index]));
                hflcal_para != nullptr)
            {
                return get_single_sigma_params(hflcal_para);
            }
            if (auto* ihflcal_para =
                    std::get_if<IHFLCALPara>(std::addressof(simplified_model_paras[index]));
                ihflcal_para != nullptr)
            {
                return get_double_sigma_params(ihflcal_para);
            }
            if (auto* bhflcal_para =
                    std::get_if<BHFLCALPara>(std::addressof(simplified_model_paras[index]));
                bhflcal_para != nullptr)
            {
                return get_double_sigma_params(bhflcal_para);
            }
            if (auto* neg_para =
                    std::get_if<NEGPara>(std::addressof(simplified_model_paras[index]));
                neg_para != nullptr)
            {
                return get_double_sigma_params(neg_para);
            }
        }

        if (params->model_type_ == ModelType::kConvolution)
        {
            if (auto* unizar_para =
                    std::get_if<UNIZARPara>(std::addressof(convolution_model_paras[index]));
                unizar_para != nullptr)
            {
                return get_single_sigma_params(unizar_para);
            }
            if (auto* huang_para =
                    std::get_if<HuangPara>(std::addressof(convolution_model_paras[index]));
                huang_para != nullptr)
            {
                return get_double_sigma_params(huang_para);
            }
            if (auto* icauchy_para =
                    std::get_if<ICauchyPara>(std::addressof(convolution_model_paras[index]));
                icauchy_para != nullptr)
            {
                const auto diffusion_coeff =
                    icauchy_para->*solar::getParaPtr<ICauchyPara, diffusion_coeff_t>();
                const auto flux_param =
                    icauchy_para->*solar::getParaPtr<ICauchyPara, flux_param_t>();
                return {diffusion_coeff, diffusion_coeff, flux_param};
            }
            if (auto* unizar2_para =
                    std::get_if<UNIZARPara2>(std::addressof(convolution_model_paras[index]));
                unizar2_para != nullptr)
            {
                return get_double_sigma_params(unizar2_para);
            }
        }
        return {};
    };

    static const auto get_convolution_area =
        [&](int index) -> std::tuple<double, double, double, double>
    {
        if (params->model_type_ == ModelType::kConvolution)
        {
            if (auto* unizar2_para =
                    std::get_if<UNIZARPara2>(std::addressof(convolution_model_paras[index]));
                unizar2_para != nullptr)
            {
                return get_convolution_area_rectangle(unizar2_para);
            }
        }
        return {};
    };

    float* d_flux_map_sum_result = nullptr;
    Int2 receiver_resolution{};

    auto* d_shadow_block_center_bias_factor = shadow_block_calculator.getResult();
    if (params->model_type_ == ModelType::kSimplified)
    {
        solar::cuda::gpu2cpuAsync(h_shadow_block_center_bias_factor_ptr,
                                  d_shadow_block_center_bias_factor, rectangle_heliostats.size(),
                                  shadow_block_calculator.getCudaStream());
    }

    if (heliostat_tmp_geometry->getType() == GeometryType::kGeometryRectangle3D)
    {
        if (params->model_type_ == ModelType::kSimplified)
        {
            if (receiver_geometry->getType() == GeometryType::kGeometryCylinder3D)
            {
                cuda::inter_op::CylinderReceiver cylinder_receiver(simulation_stream);
                cuda::inter_op::convertReceiver(std::addressof(params->receivers_[0]),
                                                std::addressof(cylinder_receiver));
                cylinder_receiver.setNormal({0, 0, -1});
                cylinder_receiver.setSurfaceIndex(0);
                cylinder_receiver.init(1.0 / receiver_pixel_size.x());
                receiver_resolution = cylinder_receiver.getResolution();
                CHECK_CUDA_ERRORS(
                    cudaMallocAsync(&d_flux_map_sum_result,
                                    sizeof(float) * receiver_resolution.x * receiver_resolution.y,
                                    simulation_stream));
                CHECK_CUDA_ERRORS(
                    cudaMemsetAsync(d_flux_map_sum_result, 0,
                                    sizeof(float) * receiver_resolution.x * receiver_resolution.y,
                                    simulation_stream));
                preprocessing_thread_pool.join();
                CHECK_CUDA_ERRORS(cudaStreamSynchronize(shadow_block_calculator.getCudaStream()));

                if (params->simulation_option_ == SimulationOption::PerHeliostat)
                {
                    // cuda::inter_op::SimplifiedMethod* previous_method = nullptr;
                    // auto start = std::chrono::system_clock::now();
                    // double obj_init_time = 0;

                    for (int i = 0; i < params->heliostats_.size(); i++)
                    {
                        // auto current_start = std::chrono::system_clock::now();
                        const auto& rectangle_heliostat = rectangle_heliostats[i];
                        const auto focus_center =
                            cylinder_receiver.getFocusCenter(rectangle_heliostat.getPosition());
                        auto current_method = cuda::inter_op::SimplifiedMethod(
                            rectangle_heliostat.getVertices().data(), cylinder_receiver.getSize(),
                            rectangle_heliostat.getNormal(), cylinder_receiver.getPosition(),
                            rectangle_heliostat.getPosition(), focus_center,
                            rectangle_heliostat.getSize(), cylinder_receiver.getResolution(),
                            SOLAR_GET_VALUE(ScenePara, DNI_, params->scene_para_),
                            d_flux_map_sum_result, simulation_stream);
                        // current_method->setCudaStream(simulation_stream);

                        const auto [sigma_x, sigma_y, peak_flux_param] = get_sigmas(i);
                        current_method.setModelParams(sigma_x, sigma_y, peak_flux_param);
                        current_method.setCenterBias({h_shadow_block_center_bias_factor[i].x,
                                                      h_shadow_block_center_bias_factor[i].y,
                                                      h_shadow_block_center_bias_factor[i].z});
                        current_method.setShadowBlockFactor(h_shadow_block_center_bias_factor[i].w);
                        // auto current_end = std::chrono::system_clock::now();
                        // auto current_duration =
                        //     std::chrono::duration_cast<std::chrono::microseconds>(current_end -
                        //                                                           current_start);
                        // obj_init_time += double(current_duration.count()) *
                        //                  std::chrono::microseconds::period::num /
                        //                  std::chrono::microseconds::period::den;
                        current_method.drawFluxMapCylinder();
                        // previous_method = current_method;
                    }
                    // not_finished = false;
                    // cudaStreamSynchronize(simulation_stream);
                    // delete previous_method;
                    // auto end = std::chrono::system_clock::now();
                    // auto duration =
                    //     std::chrono::duration_cast<std::chrono::microseconds>(end - start);
                    // std::cout << solar::format("simulation_for: {}s\n",
                    //                            double(duration.count()) *
                    //                                std::chrono::microseconds::period::num /
                    //                                std::chrono::microseconds::period::den);
                    // std::cout << solar::format("obj_init_time: {}s\n", obj_init_time);
                }
                if (params->simulation_option_ == SimulationOption::Whole)
                {
                    std::vector<Float3> flux_params(rectangle_heliostats.size());
                    Float3* d_flux_params = nullptr;
                    preprocessing_thread_pool.join();
                    SOLAR_OMP_PARALLEL_FOR_SIMD
                    for (int i = 0; i < rectangle_heliostats.size(); i++)
                    {
                        const auto [sigma_x, sigma_y, peak_flux_param] = get_sigmas(i);
                        flux_params[i] = Float3(sigma_x, sigma_y, peak_flux_param);
                    }
                    solar::cuda::cpu2gpuAsync(d_flux_params, flux_params.data(),
                                              rectangle_heliostats.size(), simulation_stream);
                    // auto start = std::chrono::system_clock::now();
                    auto current_method = cuda::inter_op::SimplifiedMethod(
                        d_image_plane_axes, d_shadow_block_center_bias_factor, d_flux_params,
                        rectangle_heliostats.size(), cylinder_receiver.getSize(),
                        cylinder_receiver.getPosition(), cylinder_receiver.getResolution(),
                        simulation_stream);
                    // cudaStreamSynchronize(current_method.getCudaStream());
                    // CHECK_CUDA_ERRORS(cudaGetLastError());
                    d_flux_map_sum_result = current_method.getFluxMap();
                    // auto end = std::chrono::system_clock::now();
                    // auto duration =
                    //     std::chrono::duration_cast<std::chrono::microseconds>(end - start);
                    // std::cout << solar::format("simulation_for: {}s\n",
                    //                            double(duration.count()) *
                    //                                std::chrono::microseconds::period::num /
                    //                                std::chrono::microseconds::period::den);
                }
                if (receiver_geometry->getType() == GeometryType::kGeometryRectangle3D)
                {
                    // cuda::inter_op::SimplifiedMethod* previous_method = nullptr;
                    cuda::inter_op::RectangleReceiver rectangle_receiver;
                    cuda::inter_op::convertReceiver(std::addressof(params->receivers_[0]),
                                                    std::addressof(rectangle_receiver));
                    rectangle_receiver.setSurfaceIndex(0);
                    rectangle_receiver.init(1.0 / receiver_pixel_size.x());
                    receiver_resolution = rectangle_receiver.getResolution();
                    CHECK_CUDA_ERRORS(cudaMallocAsync(&d_flux_map_sum_result,
                                                      sizeof(float) * receiver_resolution.x *
                                                          receiver_resolution.y,
                                                      simulation_stream));
                    CHECK_CUDA_ERRORS(cudaMemsetAsync(d_flux_map_sum_result, 0,
                                                      sizeof(float) * receiver_resolution.x *
                                                          receiver_resolution.y,
                                                      simulation_stream));
                    preprocessing_thread_pool.join();
                    CHECK_CUDA_ERRORS(
                        cudaStreamSynchronize(shadow_block_calculator.getCudaStream()));
                    for (int i = 0; i < rectangle_heliostats.size(); i++)
                    {
                        const auto& rectangle_heliostat = rectangle_heliostats[i];
                        const auto focus_center =
                            rectangle_receiver.getFocusCenter(rectangle_heliostat.getPosition());
                        const auto [sigma_x, sigma_y, peak_flux_param] = get_sigmas(i);
                        const auto receiver_normal =
                            cpu::toFloat3((cpu::toVec3(rectangle_heliostat.getPosition()) -
                                           cpu::toVec3(focus_center))
                                              .get_unit_vector());
                        auto current_method = cuda::inter_op::SimplifiedMethod(
                            rectangle_heliostat.getVertices().data(),
                            rectangle_receiver.getVertices().data(),
                            rectangle_heliostat.getNormal(), receiver_normal,
                            rectangle_heliostat.getPosition(), focus_center,
                            rectangle_heliostat.getSize(), rectangle_receiver.getResolution(),
                            SOLAR_GET_VALUE(ScenePara, DNI_, params->scene_para_),
                            simulation_stream);
                        current_method.setModelParams(sigma_x, sigma_y, peak_flux_param);
                        current_method.setCenterBias({h_shadow_block_center_bias_factor[i].x,
                                                      h_shadow_block_center_bias_factor[i].y,
                                                      h_shadow_block_center_bias_factor[i].z});
                        current_method.setShadowBlockFactor(h_shadow_block_center_bias_factor[i].w);
                        current_method.drawFluxMap();
                        // previous_method = current_method;
                    }
                    // cudaStreamSynchronize(previous_method->getCudaStream());
                }
            }

            // solar::io::csv::writeFluxMap(
            //     R"(D:\Learn\solar_software\solar_3\BHFLCAL_result_gpu_sum.txt)",
            //     d_flux_map_sum_result, receiver_resolution.y, receiver_resolution.x, 10);
        }
        if (params->model_type_ == ModelType::kConvolution)
        {
            if (receiver_geometry->getType() == GeometryType::kGeometryCylinder3D)
            {
                cuda::inter_op::CylinderReceiver cylinder_receiver(simulation_stream);
                cuda::inter_op::convertReceiver(std::addressof(params->receivers_[0]),
                                                std::addressof(cylinder_receiver));
                cylinder_receiver.setNormal({0, 0, -1});
                cylinder_receiver.setSurfaceIndex(0);
                cylinder_receiver.init(1.0 / receiver_pixel_size.x());
                receiver_resolution = cylinder_receiver.getResolution();
                // std::cout << solar::format("resolution: ({}, {})\n", receiver_resolution.x,
                // receiver_resolution.y);
                CHECK_CUDA_ERRORS(
                    cudaMallocAsync(&d_flux_map_sum_result,
                                    sizeof(float) * receiver_resolution.x * receiver_resolution.y,
                                    simulation_stream));
                CHECK_CUDA_ERRORS(
                    cudaMemsetAsync(d_flux_map_sum_result, 0,
                                    sizeof(float) * receiver_resolution.x * receiver_resolution.y,
                                    simulation_stream));
                preprocessing_thread_pool.join();
                CHECK_CUDA_ERRORS(cudaStreamSynchronize(shadow_block_calculator.getCudaStream()));
                // cuda::inter_op::ConvolutionMethod* previous_method = nullptr;

                // float* d_flux_map_sum_result2 = nullptr;
                // CHECK_CUDA_ERRORS(
                //     cudaMallocAsync(&d_flux_map_sum_result2,
                //                     sizeof(float) * receiver_resolution.x * receiver_resolution.y,
                //                     simulation_stream));
                // CHECK_CUDA_ERRORS(
                //     cudaMemsetAsync(d_flux_map_sum_result2, 0,
                //                     sizeof(float) * receiver_resolution.x * receiver_resolution.y,
                //                     simulation_stream));
                if (params->simulation_option_ == SimulationOption::PerHeliostat)
                {
                    for (int i = 0; i < params->heliostats_.size(); i++)
                    {
                        // solar::println(solar::format("{}", i));
                        // auto start = std::chrono::system_clock::now();
                        const auto& rectangle_heliostat = rectangle_heliostats[i];
                        const auto focus_center =
                            cylinder_receiver.getFocusCenter(rectangle_heliostat.getPosition());
                        const auto [sigma_x, sigma_y, peak_flux_param] = get_sigmas(i);
                        auto current_method = cuda::inter_op::ConvolutionMethod(
                            rectangle_heliostat.getVertices().data(),
                            rectangle_heliostat.getNormal(), rectangle_heliostat.getPosition(),
                            cylinder_receiver.getPosition(), rectangle_heliostat.getSize(),
                            cylinder_receiver.getSize(), focus_center, helio_row_col,
                            cylinder_receiver.getResolution(),
                            shadow_block_calculator.getMicroHeliostatShadowBlock(),
                            std::get<ConvolutionModelType>(params->detailed_model_type_),
                            d_flux_map_sum_result, simulation_stream);
                        if (use_dnn_based_model)
                        {
                            const auto [x3, x0, y1, y0] = get_convolution_area(i);
                            dnn::UNIZAR2Output model_output;
                            model_output.peak_flux_param_ = peak_flux_param;
                            model_output.sigma_x_ = sigma_x;
                            model_output.sigma_y_ = sigma_y;
                            model_output.x3_ = x3;
                            model_output.x0_ = x0;
                            model_output.y1_ = y1;
                            model_output.y0_ = y0;
                            current_method.setModelParams(model_output);
                        }
                        else
                        {
                            current_method.setModelParams(sigma_x, sigma_y, peak_flux_param);
                        }
                        current_method.setHeliostatID(i);
                        current_method.drawFluxMapCylinderFast();
                        // previous_method = current_method;

                        // auto end = std::chrono::system_clock::now();
                        // auto duration =
                        //     std::chrono::duration_cast<std::chrono::microseconds>(end - start);
                        // std::cout << solar::format("UNIZAR per heliostat: {}s\n",
                        //                            double(duration.count()) *
                        //                                std::chrono::microseconds::period::num /
                        //                                std::chrono::microseconds::period::den);
                    }
                    // for (int i = 0; i < 1; i++)
                    // {
                    //     // solar::println(solar::format("{}", i));
                    //     // auto start = std::chrono::system_clock::now();
                    //     const auto& rectangle_heliostat = rectangle_heliostats[i];
                    //     const auto focus_center =
                    //         cylinder_receiver.getFocusCenter(rectangle_heliostat.getPosition());
                    //     const auto peak_flux_param = peakFluxParam(
                    //         focus_center, rectangle_heliostat.getPosition(),
                    //         rectangle_heliostat.getNormal(), rectangle_heliostat.getSize(),
                    //         1000.0);
                    //     double sigma_x = kUndefined;
                    //     double sigma_y = kUndefined;
                    //     std::tie(sigma_x, sigma_y, std::ignore) = get_sigmas(i);
                    //     auto* current_method = new cuda::inter_op::ConvolutionMethod(
                    //         rectangle_heliostat.getVertices().data(),
                    //         rectangle_heliostat.getNormal(), rectangle_heliostat.getPosition(),
                    //         cylinder_receiver.getPosition(), rectangle_heliostat.getSize(),
                    //         cylinder_receiver.getSize(), focus_center, helio_row_col,
                    //         cylinder_receiver.getResolution(),
                    //         shadow_block_calculator.getMicroHeliostatShadowBlock(),
                    //         std::get<ConvolutionModelType>(params->detailed_model_type_),
                    //         HeliostatType::kWithoutGlass, d_flux_map_sum_result2,
                    //         simulation_stream);
                    //     current_method->setModelParams(sigma_x, sigma_y, peak_flux_param);
                    //     current_method->setHeliostatID(i);
                    //     current_method->drawFluxMapCylinderFast();
                    //     previous_method = current_method;

                    //     // auto end = std::chrono::system_clock::now();
                    //     // auto duration =
                    //     //     std::chrono::duration_cast<std::chrono::microseconds>(end -
                    //     start);
                    //     // std::cout << solar::format("UNIZAR per heliostat: {}s\n",
                    //     //                            double(duration.count()) *
                    //     //                                std::chrono::microseconds::period::num
                    //     /
                    //     // std::chrono::microseconds::period::den);
                    // }
                    // std::cout << solar::cuda::metrics::WRAD(
                    //                  d_flux_map_sum_result2, d_flux_map_sum_result,
                    //                  cylinder_receiver.getResolution(), simulation_stream)
                    //           << '\n';
                    // // CHECK_CUDA_ERRORS(cudaStreamEndCapture(simulation_stream, &graph));
                    // // CHECK_CUDA_ERRORS(cudaGraphInstantiate(&instance, graph, nullptr, nullptr,
                    // // 0)); CHECK_CUDA_ERRORS(cudaGraphLaunch(instance, simulation_stream));
                }
                if (params->simulation_option_ == SimulationOption::Whole)
                {
                    std::vector<Float3> flux_params;
                    std::vector<Float3> focus_centers(rectangle_heliostats.size());
                    Float3* d_flux_params = nullptr;
                    Float3* d_focus_centers = nullptr;
                    CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(d_focus_centers),
                                                      sizeof(Float3) * rectangle_heliostats.size(),
                                                      simulation_stream));
                    preprocessing_thread_pool.join();
                    getFocusCenters(GeometryType::kGeometryCylinder3D, cylinder_receiver.getSize(),
                                    cylinder_receiver.getPosition(), {}, d_heliostat_positions,
                                    d_focus_centers, rectangle_heliostats.size(),
                                    simulation_stream);
                    if (!use_dnn_based_model)
                    {
                        flux_params.resize(rectangle_heliostats.size());
                        SOLAR_OMP_PARALLEL_FOR_SIMD
                        for (int i = 0; i < rectangle_heliostats.size(); i++)
                        {
                            const auto [sigma_x, sigma_y, peak_flux_param] = get_sigmas(i);
                            flux_params[i] = Float3(sigma_x, sigma_y, peak_flux_param);
                        }
                        solar::cuda::cpu2gpuAsync(d_flux_params, flux_params.data(),
                                                  rectangle_heliostats.size(), simulation_stream);
                    }
                    // auto start = std::chrono::system_clock::now();
                    // cudaStreamSynchronize(simulation_stream);
                    if (!use_dnn_based_model)
                    {
                        auto current_method = cuda::inter_op::ConvolutionMethod(
                            d_image_plane_axes, d_all_heliostat_vertices,
                            shadow_block_calculator.getMicroHeliostatShadowBlock(),
                            shadow_block_calculator.getMicroHeliostatShadowBlockPerRow(),
                            d_flux_params, d_focus_centers, params->heliostats_.size(),
                            cylinder_receiver.getSize(), cylinder_receiver.getPosition(),
                            helio_row_col, cylinder_receiver.getResolution(),
                            std::get<ConvolutionModelType>(params->detailed_model_type_),
                            d_flux_map_sum_result, simulation_stream);
                    }
                    else
                    {
                        auto* const d_unizar2_outputs =
                            std::bit_cast<dnn::UNIZAR2Output*>(d_dnn_outputs);
                        auto current_method = cuda::inter_op::ConvolutionMethod(
                            d_image_plane_axes, d_all_heliostat_vertices,
                            shadow_block_calculator.getMicroHeliostatShadowBlock(),
                            shadow_block_calculator.getMicroHeliostatShadowBlockPerRow(),
                            d_unizar2_outputs, d_focus_centers, params->heliostats_.size(),
                            cylinder_receiver.getSize(), cylinder_receiver.getPosition(),
                            helio_row_col, cylinder_receiver.getResolution(),
                            std::get<ConvolutionModelType>(params->detailed_model_type_),
                            d_flux_map_sum_result, simulation_stream);
                    }
                    // cudaStreamSynchronize(simulation_stream);
                    // CHECK_CUDA_ERRORS(cudaGetLastError());
                    // auto end = std::chrono::system_clock::now();
                    // auto duration =
                    //     std::chrono::duration_cast<std::chrono::microseconds>(end - start);
                    // std::cout << solar::format("simulation_for: {}s\n",
                    //                            double(duration.count()) *
                    //                                std::chrono::microseconds::period::num /
                    //                                std::chrono::microseconds::period::den);
                }
            }
            if (receiver_geometry->getType() == GeometryType::kGeometryRectangle3D)
            {
                // cuda::inter_op::ConvolutionMethod* previous_method = nullptr;
                cuda::inter_op::RectangleReceiver rectangle_receiver;
                cuda::inter_op::convertReceiver(std::addressof(params->receivers_[0]),
                                                std::addressof(rectangle_receiver));
                rectangle_receiver.setSurfaceIndex(0);
                rectangle_receiver.init(1.0 / receiver_pixel_size.x());
                receiver_resolution = rectangle_receiver.getResolution();
                CHECK_CUDA_ERRORS(
                    cudaMallocAsync(&d_flux_map_sum_result,
                                    sizeof(float) * receiver_resolution.x * receiver_resolution.y,
                                    simulation_stream));
                CHECK_CUDA_ERRORS(
                    cudaMemsetAsync(d_flux_map_sum_result, 0,
                                    sizeof(float) * receiver_resolution.x * receiver_resolution.y,
                                    simulation_stream));
                preprocessing_thread_pool.join();
                CHECK_CUDA_ERRORS(cudaStreamSynchronize(shadow_block_calculator.getCudaStream()));
                for (int i = 0; i < rectangle_heliostats.size(); i++)
                {
                    const auto& rectangle_heliostat = rectangle_heliostats[i];
                    const auto focus_center =
                        rectangle_receiver.getFocusCenter(rectangle_heliostat.getPosition());
                    std::cout << solar::format("({}, {}, {})\n", focus_center.x, focus_center.y, focus_center.z);
                    const auto [sigma_x, sigma_y, peak_flux_param] = get_sigmas(i);
                    const auto receiver_normal = cpu::toFloat3(
                        (cpu::toVec3(rectangle_heliostat.getPosition()) - cpu::toVec3(focus_center))
                            .get_unit_vector());
                    auto current_method = cuda::inter_op::ConvolutionMethod(
                        rectangle_heliostat.getVertices().data(), rectangle_heliostat.getNormal(),
                        rectangle_heliostat.getPosition(), focus_center, helio_row_col,
                        rectangle_receiver.getResolution(), rectangle_receiver.getVertices().data(),
                        receiver_normal, shadow_block_calculator.getMicroHeliostatShadowBlock(),
                        d_flux_map_sum_result, simulation_stream);
                    if (use_dnn_based_model)
                    {
                        const auto [x3, x0, y1, y0] = get_convolution_area(i);
                        dnn::UNIZAR2Output model_output;
                        model_output.peak_flux_param_ = peak_flux_param;
                        model_output.sigma_x_ = sigma_x;
                        model_output.sigma_y_ = sigma_y;
                        model_output.x3_ = x3;
                        model_output.x0_ = x0;
                        model_output.y1_ = y1;
                        model_output.y0_ = y0;
                        current_method.setModelParams(model_output);
                    }
                    else
                    {
                        current_method.setModelParams(sigma_x, sigma_y, peak_flux_param);
                    }
                    current_method.setHeliostatID(i);
                    current_method.drawFluxMapFast();
                    // previous_method = current_method;
                }
            }
            // CHECK_CUDA_ERRORS(cudaStreamSynchronize(simulation_stream));
            // solar::io::csv::writeFluxMap(
            //     R"(D:\Learn\solar_software\solar_3\UNIZAR_result_gpu_sum.txt)",
            //     d_flux_map_sum_result, receiver_resolution.y, receiver_resolution.x, 10);
        }
    }
    grid.clear();

    if (params->simulation_option_ == SimulationOption::Whole)
    {
        if (params->model_type_ == ModelType::kConvolution)
        {
            CHECK_CUDA_ERRORS(cudaFreeAsync(d_all_heliostat_vertices, simulation_stream));
        }
    }
    CHECK_CUDA_ERRORS(cudaFreeAsync(d_heliostat_positions, simulation_stream));

    if (stream == nullptr)
    {
        cudaStreamSynchronize(simulation_stream);
        cudaStreamDestroy(simulation_stream);
    }

    return {d_flux_map_sum_result, receiver_resolution};
}

auto solar::saveFluxMap(const std::string_view& file_name, float* data, int height, int width,
                        int precision, int rows_package) -> float
{
    std::ofstream fout(file_name.data());
    std::stringstream ss;

    int address = 0;
    float max_value = 0.0F;
    // int max_value_located_row_index = 0;
    // int max_value_located_col_index = 0;
    float summation_value = 0.0F;

    for (int r = 0; r < height; ++r)
    {
        if (r % rows_package == rows_package - 1)
        {
            fout << ss.rdbuf();
            ss.clear();
        }

        for (int c = 0; c < width; ++c)
        {
            address = (height - 1 - r) * width + c;
            //            address = height * width + c;
            if (data[address] < cuda::inter_op::kEpsilon)
            {
                ss << 0;
            }
            else
            {
                max_value = std::max(max_value, data[address]);
                // if (max_value == data[address])
                // {
                //     max_value_located_row_index = r;
                //     max_value_located_col_index = c;
                // }

                summation_value += data[address];
                ss << std::fixed << std::setprecision(precision) << data[address];
            }

            if (c != width - 1)
            {
                ss << ',';
            }
            else
            {
                ss << '\n';
            }
        }
    }

    // std::cout << "\tMax value is " << max_value << ". Sum value is " << summation_value/1000<<
    // "." << std::endl; std::cout << "\tMax value is located at (" << max_value_located_row_index
    // << ", " << max_value_located_col_index << ")\n";
    fout << ss.rdbuf();
    fout.close();

    return summation_value;
}
